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Based on the conventional energy band theory, an approach is presented to describe the electronic 
structure of crystalline insulators in the presence of a finite homogeneous electric field. The expres- 
sion of polarization is derived which extends the "Berry-phase" theory of polarization to nonzero 
fields. The characteristics of the solutions are studied in details and the associations with the existing 
perturbation theories and computational schemes are discussed. 
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In recent years, much attention was paid to the the- 
ory of macroscopic polarization and the effect of homoge- 
neous electric field in crystalline solids. [3- Q Although 
the energy band theory and the density-functional the- 
ory (DFT) have provided a solid basis in understand- 
ing the equilibrium properties of crystals, the problem of 
crystalline insulators in the presence of a homogeneous 
macroscopic electric field was not fully solved. 013 The 
obstacle comes from the fact that the electric potential 
e£ ■ r is a linear term in the spatial coordinates, which 
violates the periodicity condition underlying Bloch's the- 
orem. Moreover, strictly speaking, there is no ground 
state for a crystalline system under an external field. On 
the other hand, crystalline insulators are experimentally 
stable when the applied field is not too strong, and the 
field-dependent response of the systems (such as dielec- 
tric and piezoelectric properties) is well defined in the 
measurements. This looks like a paradox. It appears 
that some works are needed to bridge the gap between 
the theory and the experiment. Some progress within 
the framework of perturbation theory has been achieved 
in this direction, |j, 0, IS nevertheless finite electric 
fields can not be handled directly in such theories due 
to the weakness of the perturbation approach. Compu- 
tational schemes applicable at finite electric fields were 
developed. H, 0, 0] In such schemes, the non-periodic 
electric potential term e£ ■ (r) in the energy functional is 
replaced by — £ ■ P where the polarization P is expressed 
in the Berry-Phase approach that holds the periodicity 
condition. They are effective in numerical computing, 
but the physical basis is not very clear. 

On the basis of a band-like electronic-structure theory 
under finite electric fields, this paper aims at providing 
basic theoretical understandings on the subject. The ex- 
pression of macroscopic polarization is derived to extend 
the "Berry-phase" theory of polarization into the cases 
of nonzero fields. It is shown that the formalism can be 
transferred into a conventional eigen problem in (r, k)- 
space. The properties of the solutions are identified, and 
the connections with the existing perturbation theories 
and computational schemes are discussed. 

The theory can be developed as what were done in 
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FIG. 1: Schematic graphics of the change of the phase E^t 
when the system evolves from the state (k, t) to (k + dk, 
t + dt). If we assume the phase keeps unchanged when the 
system transit from k to k + dk (process of dash line), the 
total variation of phase is E^dt. 



Refs. 5 and 11, or by use of the technique of the crystal 
momentum representation Here we present a simpli- 
fied derivation that keeps the insights in physics. 

Our approach is based on the conventional energy band 
theory. According to the Bloch's theorem, the eigen-state 
of an electron moving in a periodic potential V(r) can be 
expressed as 



V> k (r,t) = e" 



i-Ekt/Hgik-: 



(1) 



where «k(r) is a periodic function. When there exists an 
external electric field £, the Hamiltonian of the system 
is given as: 



H 



2m 



V 2 + V(r) +e£-r. 



(2) 



In this case, k is no longer the good quantum number 
to specify the eigen-state. In the general solid state 
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physics. pj| the response of solids to an electric field is 
described as the constant-velocity movement of the wave 
vector k, i.e., 



dk 
n— 

dt 



-e£. 



(3) 



The bands are either fully filled or empty in insulators, 
so a static equilibrium state is reached regardless of the 
movement of k. Since k varies with the time t, the dif- 
ferential of Eq. (1) with respect to t yields: 



E k - Tl— ■ r 

dt 



dt 



M k (r). 

(4) 

A key point in deriving the above equation is that we 
assume the phase E^t keeps unchanged when the system 
transits between different k states, so Vk-Ek is ignored in 
the analysis (see Fig. 1). It can be demonstrated more 
clearly as follows: if the system transits from an initial 
wave-function -0 k to a final one t/Vj then for a second 
initial wave-function ^ k exp(i(p ), the system will tran- 
sit into the final one ip k i exp(z<p ), i.e., the initial phase 
keeps unchanged. Such a process is invariant for time 
translation. In contrast, if the term iVk-Ek • dk in Fig. 1 
(dash line) is included, the symmetry of the system for 
time translation will be violated. 

Based on the above analysis and the Schrodinger's 
equation ihdip k /dt = Htjj^, we get the main equation 
of crystalline insulators under electric field as: 



(H + ie£ ■ V k )uk = E k u k {r), 



(5) 



where 
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+ V(r). (6) 



It can be seen that the introduction of an electric field 
causes interactions between different k states, and Uk(r) 
can be no longer solved independently. However, the pe- 
riodicity condition in r-space can be safely applied upon 
Eq. (5), while Eq. (2) is not periodic, thus it overcomes 
the difficulties mentioned above. Essentially, Eq. (5) is 
equivalent to what is obtained in Refs. 5 and 11, but it 
is more convenient for the following analyzing. 

To one's surprise, the calculation of polarization in 
electronic structure theory remained a challenge until the 
Berry-phase theory of polarization was developed in the 
early 1990s. Q, 0, El Although such polarization theory 
has been successfully applied in many cases, its validity 
under nonzero electric field was not proven in any publi- 
cations. With Eq. (5), the expression of polarization at 
finite fields can be derived easilyQ . 

Consider an infinitesimal variation of electric field, 6£, 
in the system. The linear terms of Eq. (5) with respect 
to S£ yields 



{H + ie£ ■ V k 



S£ 



zeV k u k = E k 



5E V 



S£ S£ 



u k . (7) 



Multiplying the above equation by « k and integrating it, 
one gets 



SE k 
S£ 



dk = I ie(itk| Vk|wk}dk 



(8) 



if the periodic condition is adopted in k-space as 

w k+K (r) = e-^ r u k (r), (9) 

where K is any reciprocal lattice vector. So the macro- 
scopic polarization is achieved as 



(5¥F / E * dk 



S£ 



(2nf 



W|V k |it k )dk, 



(10) 



which is consistent with the result given by the Berry- 
phase theory. 0, 0] Although the original derivation of 
the polarization's formula in the Berry-phase approach 
assumed £ = 0, the above analysis showed that the same 
formula keeps valid for f ^ 0. 

Gauge transformation plays an important role in 
Eq. (5). It is noted that Eq. (5) is invariant under the 
following transformation: 



u k = Mkexp [i(3(k)] , 
E k = E k - e£ ■ V k /3(k), 



(11) 



where /?(k) is an arbitrary function. Such phase arbi- 
trariness can not be removed even if the periodic bound- 
ary condition [Eq. (9)] is adopted. This reveals that the 
polarization in Eq. (10) is a multivalued quantity. 

Gauge transformation will not change the physical 
phase of the system. Within an infinitesimal time in- 
terval dt, the extra phase of itk cased by gauge trans- 
formation is 8ip 1 — Vk/?(k) • dk while that of e ^ lEkt ^ h is 
Sip 2 = e£-Vk[3(k)dt/h. With (3) we arrive 5(p 1 +5ip 2 — 0. 
So the gauge transformation just moves a part of the 
phase from Uk to e~ lEkt l h while keep the total phase 
of u-k.e~ lEkt l h unchanged. It does not cause any mea- 
surable effects. Furthermore, such an arbitrariness can 
be utilized to simplify some analyses. By choosing /3(k) 
appropriately, E^ can be transformed into a constant in- 
dependent of k, i.e., -Ek = E, and thus Eq. (5) becomes 
a conventional eigen-problem in (r, k)-space: 



{H + ie£ ■ V k )itk = Eu k (r), 



(12) 



whose properties can be well understood based on the 
existing knowledge of eigen problem. Some properties 
of the solution are listed in the following: (1) if there 
are N bands for £ = 0, then Eq. (12) has N physically 
independent solutions for £ ^ 0, each of which is as- 
sociated with a group of physically equivalent solutions 
as Uk = Uk exp (ik ■ R) and E — E — e£ ■ R where R 
is any lattice site; (2) Vk(uk|tik) = so that Uk can 
be normalized as (ukl^k) = 1, which is not self-evident 
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FIG. 2: (a) Band energy E, (b) polarization P and (c) usual 
energy with polarization's contribution excluded, Eo — E + 
£ ■ P, of three bands for the tigh-binding model (see Eq. (14)) 
as functions of electric field. Labels I, m and n in Fig. 2(c) are 
used to specify the index n of peaks between different bands 
(see Eq. (16)). 



for £ 7^ because different k states are correlated by 
Eqs (5) or (12); (3) the wavefunctions coming from dif- 
ferent bands and with the same k value are orthogonal, 
i.e., (Wkl u k) = ' wnere a an d (3 denote different bands. 

Eq. (12) is equivalent to the variational analysis on the 
energy functional: 



(E) 



(2 



^Tg J (u k \H \u k )dk-S-P 



(13) 



with P in Eq. (10), which is the basis of the perturbation 
theories and computational schemes in references 0, H, 

EJ. 

To learn more characteristics of the solution, we con- 
duct calculations on a ID tight-binding Hamiltonian with 
three bands: 81 



H 



3 C j C 3 



+ t 



Cj+l + c j + l c j 



(14) 



where the site energy is given as e^m+k = Acos(a — 
27rfc/3). a, t and A are three parameters of the Hamilto- 
nian. The position operator is set to be x = X)j(j/3) c j c j 
since there are three atoms in each unit cell. 

In Fig. 2, the properties of three bands are reported as 
functions of electric field for the parameter e = t = A = 1 
and a = 0. The most extraordinary feature at finite elec- 
tric fields is that there are many peaks in the curves of 
polarization P and the usual energy with polarization's 
contribution excluded, E (E — {H ) = E + £ ■ P). 
Correspondingly, the curves of band energy E are dis- 
continuous with many gaps. 

The peaks are caused by the overlap of the energy of 
different bands. According to Eq. (12), the dielectric 
susceptibility of the band a is given as: 



2e 2 



' (2irfe 



(T/lVkK^t/lVklu") 
E0 - E a ' 



Note that the band energy is a multivalued quantity with 
a term —ne£a where a is the crystal lattice in the current 
one-dimensional case. So if the electric field satisfies the 
following relation: 



1 E 13 - E° 

£ = 

n ea 



(n = ±l,±2,...), 



(16) 



the dielectric susceptibility diverges and its sign changes 
when the electric field goes across such points, which re- 
sults in the peaks of the polarization. 

The revealed properties of polarization at finite fields 
here are helpful in understanding the convergency of the 
perturbation theories. The electric filed associated with 
the peaks of polarization is a kind of "singular point" . 
From Eq. (16) it is known that there are more and more 
such singular points when the electric filed is getting 
smaller. Especially, the zero field £ = is a non-isolated 
singular point of which there are always other singular 
points within any small neighborhoods. Therefore, if 
the properties of the system are expanded in the form of 
power series in the electric field, the series are not con- 
vergent within any finite neighborhood of radius centered 
about £ = 0. That is the reason why the perturbation 
theories in continuum formulation are valid only for in- 
finitesimal fields. 0,111 

When the band energy E has a minimal value, the so- 
lution of ground state can be obtained by minimizing the 
energy functional Eq. (13). Unfortunately, the minimal 
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FIG. 3: (a) The average polarization and (b) dielectric sus- 
ceptibility of the three-band model. The peaks in Fig. 2 have 
been smoothed out to provide results comparable with exper- 
iments. 



value of E does not exist in continuum formulation due 
to the term —ne£a. A technique to remove such diffi- 
culty is to perform a discretization on the k-space, i.e., 
to restrict the maximal value of the number n, which has 
been adopted in numerical schemes.^. I id] However, when 
the electric field is large enough to shift the higher band 
to below the lower band, i.e., £ > (E? - E a )/(N k ea) 
where is the discretization number in k-space, the 
minimizing procedure will fail in obtaining the ground 
state. That is the phenomena of collapse appeared in 
the numerical schemes. It is clear that such phe- 



nomena are associated with the breaking of the numerical 
schemes, but not with the breaking of the system such 
as Zener tunneling^}. 

Although there are infinite peaks in the calculated 
curve of polarization around £ = 0, the properties of 
the system is smooth enough in experiments. So it is 
necessary to have a short discussion on the width of the 
peaks. After a simple analysis of degenerate perturba- 
tion around the peak, it is shown that the width of the 
peaks, AS, is approximately estimated as: 



a 3 e x 



1 1/2 



2« - ES) 



(17) 



where Eq and Eq are the average band energy when 
there is no electric field. For a numerical estimate, we 
set x = 10, a = 2k and E$ - Eg = 2eV. Then we get 
[a 3 x/(2(£^ - Eg))] 1 / 2 = 3.3 x 1CT 11 (in SI unit). It is 
obvious that such narrow peaks at small electric fields 
are too difficult to be observed experimentally. There- 
fore, the properties obtained in experiments reflect an 
average effect where the peaks are smoothed out. In 
Fig. 3 we depict the curves of the average polarization 
and dielectric susceptibility after the peaks are removed. 
Nonlinear effect can be observed. The dielectric suscep- 
tibility is non-monotonic, which means some high-order 
susceptibility terms are negative. 

It is noted that some issues discussed in this paper, 
such as the multivalued property of the polarization and 
the absence of ground state at finite electric fields, are 
caused by the periodic boundary conditions imposed on 
the insulators. If we consider a finite piece of a solid 
and surface effects, such problems will disappear. That 
is what should be kept in mind when one considers the 
effect of electric fields on insulators. 
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APPENDIX A: RIGOROUS DERIVATION OF 
POLARIZATION 

In crystalline insulators, the sole value of polarization 
is meaningless since its value depends on the state of the 
surface. The meaningful one, which can be compared 
with experiments, is the change of polarization induced 
by external factors such as electric field or stress. So a 
rigorous derivation of polarization is based on the anal- 
ysis on the variation of the system with respect to the 
external electric field. 

Denote the expected value of H as 



E = (u\H \u). 



(Al) 



It is a physical quantity that is well defined in crystalline 
insulators. Then the physical energy of the system can 
be written as: 



(E)=Eo-P- £, 



(A2) 



where P is the physical polarization whose expression is 
unknown at this point. For a ground (stable) state, the 
energy functional reaches its minimum, i.e., 



6(E) = 



(A3) 



for any infinitesimal variation of wavefunctions. Now, 
considering the variation of wavefunctions in a virtual 
process £ — > £ + S£, from Eqs. (A2) and (A3) we have 



So 



dP dE 
d£ ' d£ ' 

It relates P to the well-defined quantity E . 
Now define the Berry-Phase polarization as: 



(A4) 



(A5) 



le 



(2tt) 3 



(uk|V k |uk)dk:. 



(A6) 



It is not a physical quantity since it changes with gauge 
transformation. From Eq. (5), we have 



£ ■ Pr — En — 



1 



(2tt) 3 



£ k dk. 



(A7) 



Conducting differential on the above equation with re- 
spect to £ , we have 



£■ 



dPv 



dE n 



E k dk 



d£ d£ (2tt) 3 d£ 
Combined with Eqs. (10), (A5) and (A6), we reach 
dPn dP 



(A8) 



(A9) 



d£ d£ 

So, the polarization defined in Berry-Phase expression, 
Pb , is equal to the physical polarization P except a con- 
stant. 



APPENDIX B: NORMALIZATION AND 
ORTHOGONALITY OF u k 

Under finite electric field, there are correlations be- 
tween wavefunctions Uk with different k values. However, 
it can be shown that Uk are normalized and orthogonal 
as the case of £ = 0. 

Compared with the Schrodingcr's equation 
ihd^ k /dt — Hip w , it is recognized that the role of 
k in Eq. (5) is similar to that of t. A consequence of 
this similarity is the conservation of possibility in the 
k-space: 



Vk(uk|«k) = 0, 



(Bl) 



which is a physical requirement for the evolution equa- 
tion. So the wavefunction can be normalized as 
(uk\uk) = 1 or (itk|«k) = a 3 - The difference be- 
tween k and t exists in the boundary condition. In 
the Schrodingcr's equation, no boundary condition is re- 
quired in i-space, while in Eq. (5) k and k + K label the 
same state (with a possible phase difference) , so that the 
solutions to Eq. (5) are discrete (band structure). 

In the conventional band theory without electric field, 
the wave-functions coming from different bands are or- 
thogonal, i.e., 



<«£!«£> = o, 



(B2) 



where a and (3 denote different bands. Since Eq. (12) is 
a conventional cigen-problem, the orthogonality gives 



E 



<«£!«£> = o 



(B3) 



at finite fields. In the first glance, Eq. (B2) is no longer 
valid. Nevertheless, to one's surprise, Eq. (B2) keeps 
valid in this case. To prove this point, we transform u^, 
into 



= exp(zk • R), 



(B4) 



(R can be any lattice vector), which is also an cigen wave- 
function of Eq. (12), thus 

E«l"k> = 5>kl«£> exp(zk ■ R) = 0. (B5) 



By applying the knowledge of Fourier's transformation 
on Eq. (18), we arrive at Eq. (B2) with electric fields. 



APPENDIX C: THE WIDTH OF THE PEAK OF 
POLARIZATION 

The width of the peak can be evaluated by considering 
the degenerate perturbation analysis on two bands u a 
and vP . Expanding the wavefunction on the bands a 
and P of £ = 0, the Hamiltonian becomes: 



H = 
where 



Eg + £- a 3 P aa £ ■ a 3 P Q/3 
£-a 3 P? a E^+£-a 3 P^ 



(CI) 



(2tt) £ 



«olVkK }dk (^i/ = a,/3). (C2) 



Eq. (CI) describes the transition between bands a and 
/?, where the transition range (width of the peak) is: 



A£ P Q/3 

£ pace _ p/3/3 ' 

The degenerate condition requires that 



Eg + £ ■ a 3 P aa = E% + £- a 3 pW, 



which gives 



p/3/3 



~ E 
£a 3 ■ 



From Eq. (15), P Q/3 is approximated as 

_ 2a 3 |P^| 2 
X ~ e {El-E§) 

Combining Eqs. (C3,C5,C6), we arrive 



(C3) 



(C4) 



(C5) 



(C6) 



A£ 



= £ 



a 3 e x 



2« - Eg) 



1/2 



(C7) 



APPENDIX D: ITERATION SOLVING OF THE 
EQUATION 

In the above main text, it is shown that our theory pro- 
vides the basis for existing computational schemes|9ill0|. 
Based on the theory, other computational schemes may 
also be developed. The calculation in the main text 
is conducted by solving the conventional eigen problem 
Eq. (12) with subroutine ZHEEVX in the program packet 
LAPACK, which is applicable for any finite electric fields 
but requires more computer capacity. Here, we demon- 
strate a second scheme to solve Eq. (5). 
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FIG. 4: Polarization of three bands given by computational 
scheme of iteration Eq. (Dl) (red), compared with what is 
given by the eigen-problem subroutine of LAPACK (blue), 
(a) Polarization of the first band, (b) Polarization of the 
second and the third bands. 



The computational scheme is based on the following 
iteration: 



(H 



o^ie£-V k \u£- 1 \r)) 
EW)u£\r)=-ie£-V k ut 1) (r) 



(Dl) 



k-space is discretized into 99 points (N^ — 99). Iteration 
stops when the discrepancy is smaller than the conver- 
gence criterion: \u k — *'| < 1CP 4 . A computation 
example on the three-band model is shown in Fig 4. It 
can be seen that the convergence is obtained at £ < 0.2 
for the first band (Fig. 4(a)), which exceeds the critical 
field given by Ref. 9 (E gap /eaNk = 0.012) by about one 
order. 



APPENDIX E: DIELECTRIC SUSCEPTIBILITY 
AT ZERO FIELD: A COMPARISON 

The analytic solution for the static susceptibility at 
zero field can be computed from the Kubo formula. In 
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the following we show that the susceptibility derived from 
the theory in this paper is consistent with what is given 
by Kubo formula. 

According to the reference (H. Ehrenreich and M. H. 
Cohen, Phys. Rev 115, 786), the analytical expression 
of the static susceptibility at zero field is given as: 



X 



AE7 h{E k i)f$M'l)- 2 d% (El) 

U17T I 
1,1' J 



where fo(Eki) is the Fermi distribution and 

2 



J i' i 



VI 



1 



?iu>i'im 



Pin = ^J u* kl ,P»u kl d 3 x, 
P>* = -ihV ■ fi. 
At zero temperature (T = 0), 

occupied 2fc , , n i > ,2 

2e h s-^ , f \{UM<\H ■ V\Ukl)\ , 3 



(E2) 

(E3) 
(E4) 



X 



' m 2 7r 2 ^— ' 



(Tfc. (E5) 



For the current theory, with an infinitesimal variation 
<5£, the linear term of Eq. (5) at £ = is: 

Expand %i as ^ = au'Uw, thus 



Integral the above equation with (ttki' | yielding 

ie 



<5£, 



fUjJl'l 

So the dielectric susceptibility is 



Wkc|V k Wk;)- 



(E7) 
(E8) 



_ 2e 2 p„V^ / {"ki ' I Vk I "kl ) * ("ki ' I Vk I "ki ) ,j3 K 

- WFn^i^i' J «,„ 01 k - 

(E9) 

Eq. (E9) can be further simplified. Multipling Vk with 
Eq. (5) gives: 



(E10) 



= ^k/VkWki + MkiVk-Eki- 
Multipling with (uw | yields 

(ukJ'|V k |uk/) = — •- — -. (Ell) 

So the dielectric susceptibility is 

(E12) 



Considering the difference of the units used (SI vs. 
Gauss), Eqs. (E12) and (E5) are consistent. 

APPENDIX F: A SHORT NOTE ON ZENER 
TUNNELING 



Essentially, the Zener tunneling is a kinetic effect in 
which the transition rate between different states should 
be calculated. Therefore, Zener tunneling can not be 
directly solved in the framework of this paper or the 
existing computational schemes 0, ITol | or perturbation 
theories |5| whose purpose is to obtain the eigen states 
of the system. Eigen states are unchanged with time and 
no transition will occur if no extra term is introduced 
into the Hamiltonian. 

Zener tunneling is usually analyzed by semi-classical 
method. Here we illustrate some properties of motion by 
the theory of this paper. 

Firstly, the velocity of the electron is given as: 



(Fl) 



Multipling Vk with Eq. (5) gives: 

H V k u k + ie£ ■ V k V k Uk - «f£Vti k + 
= -EkVkUk + WkV k £ k . 

Multipling with \uu) yields 



ih 2 ?i 2 k 

V k (M k |ie£:-Vk|Mk) (w k |V|u k ) + = V k £ k - (F3) 

m m 



(F2) 



So the velocity is 



where 



£° = (u k |£r |«k>. 



(F4) 



(F5) 



Secondly, we consider the motion of the wave packet. 
Consider a wave packet (ko,i"o) at to. Then at time t, 
the k will evolve into k = k — e£(t — t )/h, and 



r - ro = / v,II = j VkE^dt = - E * f™ . (F6) 



e£ 



to 



to 



So the at finite fields should be used to replace the 
quantity at zero field to analyze the Zener tunneling. 



